// Author: Mark Richardson
// Date: 05/02/2022
// Multi-rater item response model for latent traits with naive prior on the latent traits

// Revision date: 10/10/2022
// Revision purpose: Delete extraneous dept_group vector

// Revision date: 10/12/2022
// Revision purpose: Add pid and appointee variables to the model

data {
  
  int<lower = 0> N; // number of ratings
  
  int<lower = 0> A; // number of agencies
  
  int<lower = 0> R; // number of raters
  
  int<lower = 0> D; // number of departments
  
  int<lower = 1, upper = A> aa[N]; // agency index
  
  int<lower = 1, upper = R> rr[N]; // rater index
  
  int<lower = 1, upper = D> dd[R]; // department index
  
  vector[N] y; // vector of ratings
  
  real<lower = 0> sd_naive;
  
  vector<lower = 0, upper = 1>[N] appointee; // Appointee indicator
  
  vector<lower = 0, upper = 1>[N] rep; // Republican indicator

  vector<lower = 0, upper = 1>[N] dem; // Democrat indicator
  
}


parameters {
  
  vector[D] mu_alpha; // dept mean for rater intercepts
  
  vector<lower = 0>[D] mu_beta; // dept mean for rater slopes
  
  vector<lower = 0>[D] sigma_alpha; // dept scale for rater intercepts
  
  vector<lower = 0>[D] sigma_beta; // dept scale for rater slopes
  
  vector[A] x; // latent performance

  real<lower = 0> sigma; // scale of latent performance
  
  vector[R] alpha_tilde;
  
  vector<lower = 0>[R] beta_tilde;
  
  vector[3] gamma; // indicator coefficients
  
}

transformed parameters {
  
  vector[R] alpha; // rater intercept parameters
  
  vector[R] beta; // rater slope parameters
  
  for (r in 1:R) {
    
    alpha[r] = mu_alpha[dd[r]] + sigma_alpha[dd[r]] * alpha_tilde[r];
    
    beta[r] = mu_beta[dd[r]] + sigma_beta[dd[r]] * beta_tilde[r];
    
  }
}


model {
  
  x ~ normal(0, sd_naive); // Informed prior on latent trait
  
  mu_alpha ~ normal(0, 5);
  
  mu_beta ~ normal(0, 5);
  
  alpha_tilde ~ std_normal();
  
  beta_tilde ~ std_normal();
  
  sigma ~ normal(0, 5);
  
  sigma_alpha ~ normal(0, 5);
  
  sigma_beta ~ normal(0, 5);
  
  gamma ~ normal(0, 5);
  
  for (n in 1:N) {
    
    y[n] ~ normal(alpha[rr[n]] + beta[rr[n]] * x[aa[n]] + gamma[1] * appointee[n] + gamma[2] * rep[n] + gamma[3] * dem[n], sigma);

  }
}

